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We have recently presented a new approach for numerical relativity simulations in spherical polar coordinates, 
both for vacuum and for relativistic hydrodynamics. Our approach is based on a reference-metric formulation of 
the BSSN equations, a factoring of all tensor components, as well as a partially implicit Runge-Kutta method, 
and does not rely on a regularization of the equations, nor does it make any assumptions about the symmetry 
across the origin. In order to demonstrate this feature we present here several off-centered simulations, including 
simulations of single black holes and neutron stars whose center is placed away from the origin of the coordinate 
system, as well as the asymmetric head-on collision of two black holes. We also revisit our implementation of 
relativistic hydrodynamics and demonstrate that a reference-metric formulation of hydrodynamics together with 
a factoring of all tensor components avoids problems related to the coordinate singularities at the origin and on 
the axes. As a particularly demanding test we present results for a shock wave propagating through the origin 
of the spherical polar coordinate system. 
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I. INTRODUCTION 


Numerical relativity simulations of black holes and other 
compact objects have made remarkable progress in recent 
years. In particular, simulations of the complete binary black 
hole coalescence, from inspiral through merger to the quasi¬ 
normal ring-down of the merger remnant, became possible 
with the calculations of mm- Since then, a number of differ¬ 
ent groups have assembled accurate gravitational waveforms 
emitted in these mergers (see, e.g., the compilation by the 
NINJA collaboration, ID), and have explored astrophysical 
consequences of these mergers, including black hole recoil 
(e.g. @0) and spin flip (e.g. (8)). 

Many current numerical relativity codes (in three spatial di¬ 
mensions) share several features: they adopt the Baumgarte- 
Shapiro-Shibata-Nakamura (BSSN) formulation of Einstein’s 
equations mm, use finite-difference methods in Cartesian 
coordinates, and adopt moving puncture coordinates, i.e. a 
combination of 1+log slicing lfl2l and the “Gamma-driver” 
condition urn (a notable exception is the SpEC code; see, e.g., 
CD-) While Cartesian coordinates are well suited for many 
applications, in particular simulations of binaries, spherical 
polar coordinates have some desirable properties whenever 
the object under consideration is close to spherical or axial 
symmetry. Specific examples include gravitational collapse, 
supernova explosions, and accretion disks. 

We have therefore developed and implemented a new ap¬ 
proach that applies in spherical polar coordinates the numer¬ 
ical methods that have previously proven to be extremely 
successful in Cartesian coordinates ca As we will review 
in more detail in Section [II| this approach relies on three 
key ingredients: a reference-metric formulation of the BSSN 
equations [16, 17 ], factoring out appropriate geometrical fac¬ 
tors from tensor components, and using a ’’partially implicit” 
Runge-Kutta (PIRK) method |[T8l420ll . The resulting equa¬ 
tions are still singular at the origin of the coordinate system 
and on the polar axis, but all singular terms can be handled an¬ 
alytically, and the PIRK method is stable even in the presence 
of these singular terms. Our approach therefore does not rely 


on a regularization of the equations, and can be used even in 
the absence of spherical or axi-symmetry. In [ 211 we applied 
a reference-metric approach to the formulation of relativistic 
hydrodynamics, and implemented the resulting equations to 
perform what we believe are the first self-consistent and stable 
simulations of general relativistic hydrodynamics in dynami¬ 
cal spacetimes in spherical polar coordinates without the need 
for regularization or symmetry assumptions. 


While we have previously performed and presented a num¬ 
ber of different tests, the vast majority of those tests featured a 
symmetry about the origin, which raises the question whether 
the stability of our methods hinges on that symmetry. In 
this paper we address this issue by presenting a number of 
new tests for configurations that do not satisfy that symmetry. 
Specifically we will consider an off-centered Schwarzschild 
black hole (Section |III A 1|), an asymmetric head-on collision 
of two black holes in (Section |III A 2| ), and an off-centered 
neutron star (Section [ill B 2| ). 


The other purpose of this paper is to revisit our implemen¬ 
tation of hydrodynamics. As we will discuss in more detail in 
Section [nB| we consider here a modification to the implemen¬ 
tation that we presented in ED- As a key test we show in Sec¬ 
tion |IIIB 1| results for a shock wave that propagates through 
the origin of the coordinate system. 


We would like to emphasize that the purpose of the simula¬ 
tions presented in this paper is purely to demonstrate a point 
of principle. Placing spherical objects away from the origin 
of a spherical coordinate system clearly defeats the purpose 
of such coordinate systems from a computational perspective 
- but it does provide extremely powerful tests of the proper¬ 
ties of our computational methods. Moreover, applications for 
which we expect spherical polar coordinates to be useful, for 
example supernova collapse or accretion onto a black hole, 
may involve processes in which asymmetries move the center 
of the central object away from the origin of the coordinate 
system - it is therefore important to calibrate the performance 
of the numerical methods for such off-center configurations. 
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II. NUMERICAL RELATIVITY IN SPHERICAL POLAR 
COORDINATES 

A. Einstein’s Field Equations 


We refer the reader to m for a detailed discussion of 
our implementation of Einstein’s field equations, based on 
the BSSN formulation (9 - El |, in spherical polar coordinates. 
Here we provide a brief discussion of the main components, 
namely a reference-metric formulation of the BSSN equations 
(Section IIA 1| ), a factoring of tensor components (Section 
, and a partially-implicit Runge-Kutta scheme (Section 
. We will also include brief Sections on the gauge con¬ 


IIA 2 


II A3 


ditions used in this paper (Section II A4| ) as well as the nu¬ 
merical implementation (Section IIA 5). 


1. A Reference-Metric Formulation of the BSSN Equations 

In a 3+1 decomposition of spacetime the spacetime metric 
is written as 

ds 2 = g ab dx a dx h = —a 2 dt 2 +7 ij(dx l + /3 l dt)(dx J + ftdt), 

( 1 ) 

where a is the lapse function, /3 l the shift vector, and 

lab = gab + n a n b ( 2 ) 

the spatial metric (see, e.g., l ElGHES for textbook introduc¬ 
tions). In terms of the lapse and the shift, the normal vector 
n a on each spatial slice can be written as 

n a = (— a, 0, 0, 0) or n a = (1/a, — fd l /a). (3) 

Here and in the following indices a, b ... run over spacetime 
indices, while i, j .. . run over spatial indices only; we also use 
geometrized units with c = G = 1 throughout this paper. 

The BSSN formulation (9--1L] of Einstein’s equations fur¬ 
ther adopts a conformal rescaling of the spatial metric, 


7 ij = (4) 

where is the conformally related metric, e$ the conformal 
factor, and where we refer to <fi as the “conformal exponent”. 
This decomposition is not unique, as it allows for different 
choices of the determinant 7 of the conformally related metric, 
which then result in different values of 

e 4 ^ = ( 7 / 7 ) 1/3 . (5) 

The original BSSN formulation was based on the choice 7 = 
1, which simplifies several expressions. This choice, which 
also results in the appearance of tensors with non-zero weight, 
is appropriate in Cartesian coordinates, but not in curvilinear 
coordinates. In spherical polar coordinates one might work 
around this problem by choosing 7 = r 4 sin 2 0 instead, but a 
more elegant and powerful approach is to adopt a reference- 
metric formulation (see |[T6lfT7ll : see also 12411251). 

In a reference-metric formulation we introduce a new refer¬ 
ence metric 7 ij, together with its associated connection Y l - k . 


Strictly speaking, only a reference connection is needed for 
this formalism, but for ease of presentation we assume that 
this reference connection is associated with a reference met¬ 
ric. In principle, the reference metric could be any metric, 
but the formalism is most useful for our purposes when this 
reference metric is chosen to be the flat metric in whatever 
coordinate system is used - in our case in spherical polar co¬ 
ordinates. We then define 

Ar % j k = f) k - f) k , (6) 

where Y L - k are the connection coefficients associated with the 
conformally related metric 7 ^. As the difference between two 
connections, the coefficients ATL transform as tensors, un¬ 
like the connections themselves. We compute the coefficients 
AI 7 f rom 

Ar 7 = \ f l (pjiki + Vkiji - Viljk), (7) 

where V denotes the covariant derivative associated with the 
reference metric 7 ^. We also define the conformal connection 
functions as 

v = 7 ifcA 7 fe , (8) 

but treat these as new independent variables in the equations. 

In order to specify the conformal factor e^ we specify the 
time evolution of the determinant of the conformally related 
metric, 


da = 0 , (9) 

which Brown (H calls the “Lagrangian” choice. 

Using the above expressions, the BSSN equations for 7 ^, 0 
and other curvature quantities can be expressed independently 
of any particular choice for 7 ; see eqs. ( 21 ) in m or eqs. (9) 
in na. Moreover, many of the differential operators can now 
be expressed in terms of V. Choosing 7 ^ to be the flat metric 
in spherical polar coordinates, 

(10 0 \ 

lij = I 0 r 2 0 , ( 10 ) 

\0 0 r 2 sin 2 9 ) 

we can then express these differential operators analytically in 
terms of the spherical polar connection coefficients Y l - k . 


2. Factoring of Tensor Components 

Differential equations, when expressed in spherical polar 
coordinates, often feature singular terms at the origin or on the 
axis, where r or sin 6 vanish. The advantage of the reference- 
metric formulation of Section II! A ll is that it allows us to ex¬ 
press these singular terms analytically. However, if the vari¬ 
ables in the differential equation are tensors, then the tensor 
components may also become singular at the origin or on the 
axis. In order to treat these singular terms analytically as well, 
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we factor out appropriate powers of the geometrical factors r 
and sin 9 from tensor components. 

We write the conformal connection functions ([5]), for exam¬ 
ple, as 

A' : = f A e /r | , ( 11 ) 

y A^/(rsin#) J 

and adopt the coefficients A r , X° and A^ , which remain reg¬ 
ular in regular spacetimes, as our dynamical variables. Co¬ 
variant derivatives of A\ for example, can then be expressed 
in terms of the new variables A 2 and their derivatives. As a 
concrete example we compute 

V v K e = d v (\ e /r) + AT^ = 1 (d v \ e - cos 9X V ) , (12) 

where we have used Y e = — sin 6 cos 6. A complete list of 
all these derivatives is given in eq. (26) of urn As advertised, 
the singular behavior in the tensor components can now be 
treated analytically. 

We similarly express the conformally related metric as 

Ifij = Ifij H“ Gji 5 (13) 

where the corrections do not need to be small, and then 
write 



j hqr-y 

rh r e 

r sin 6h rif 

€ ij ~ 

rh r6 

r 2 hee 

r 2 sin#/^ 


\ r sin 0h r(p 

r 2 sin Ohoy 

r 2 sin 2 0h vv> 


Similar to our example above, the derivatives 'Dtfjk = 'Di^jk 
can then be written in terms of the variables hij - see eq. (25) 
in ca for a complete list. All other tensorial quantities are 
treated in a similar way. 

With the help of these rescalings, all variables remain fi¬ 
nite for regular spacetimes even in spherical polar coordi¬ 
nates. The equations do feature singular terms, but these sin¬ 
gular terms are treated completely analytically. We do not 
attempt to regularize the equations; instead we adopt a nu¬ 
merical method that can handle these singular but analytical 
terms. 


3. Partially Implicit Runge-Kutta 

The “partially implicit Runge-Kutta” (PIRK) method was 
introduced in fl 8 l for the BSSN equations in spherical sym¬ 
metry (see also fT9ll20lD . In particular, it was demonstrated 
that the PIRK method can handle the singular terms that ap¬ 
pear in spherical polar coordinates as long as they are treated 
analytically. We refer to the above references, as well as m, 
for a more detailed discussion of the PIRK method; here we 
illustrate the approach for a simple wave equation 

- <9 2 $ + V 2 $ = 0. (15) 

We bring this equation into a form that mimics that of the 
BSSN equations by introducing a new variable k, = — d t &. 


Also assuming spherical symmetry we then rewrite eq. © 
as a pair of first-order-in-time equations 

d t & = —k (16a) 

d t K = — — (2/r)9 r T>. (16b) 

We now recognize that the variable in the singular term, i.e. <f> 
in the term (2 /r)d r $, is evolved with an equation that does 
not feature any singular terms (eq. |16a| ). The idea is then to 
evolve eq. ( |16a| ) explicitly, and use the updated values of <I> 
when evaluating the singular terms in eq. ( |16b| ). In a fully 
implicit scheme all terms on the right-hand side of the equa¬ 
tions would be evaluated using values on the new time level, 
while in the PIRK scheme only part of the variables are eval¬ 
uated on the new time levels - namely those that appear in 
singular terms, which are also those that are evolved with a 
regular equation. The effect of this is quite dramatic: while a 
simple explicit finite-difference evolution of eqs. © quickly 
becomes unstable, this PIRK method can handle the singular 
term without problems. The advantage of PIRK over a fully 
implicit scheme is that it does not require inversion of any ma¬ 
trices; in fact, the computational cost of PIRK is very similar 
to that of fully explicit methods. In a further similarity with 
fully explicit methods, PIRK is stable only as long as the time 
step is limited by a Courant condition (see eq. ( |22| ) below). 

It turns out that the BSSN equations have a structure simi¬ 
lar to that of ( [T 6 ] ), in particular, all variables in singular terms 
obey regular equations themselves. We can therefore apply 
the PIRK method as described above (see na, including Ap¬ 
pendix B, for details.) 

4. Gauge Conditions 

We adopt different versions of “moving-puncture” coordi¬ 
nates in this paper, i.e. a combination of 1 +log slicing and 
the “Gamma-driver” condition. Specifically, we adopt both a 
“non-advective” version 

d t a = -2 olK (17) 

and an “advective” version 

d t OL - P%a = -2 aK (18) 

of 1 +log slicing [ 12 ] as a condition for the lapse function. 
Here K is the trace of the extrinsic curvature 

Kij = (19) 

where D is the covariant derivative associated with the spatial 
metric 7 ^ . We note that for stationary solutions, for which 
d t a = 0 , the non-advective condition © is consistent with 
maximal slicing K = 0. 

We also use different conditions for the shift vector /3\ The 
simplest choice is /3 l = 0, but we also use an “non-advective” 
version of the Gamma-driver condition da 

dtF = 

d t B i = n s d t A\ 


( 20 a) 

( 20 b) 
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where B l is an auxiliary vector, as well as an “advective” ver¬ 
sion of a related condition 

dtP ~ P j djp = /x 5 A < (21) 

(see, e.g., [26]). We use ps = 3/4 in both conditions. 

5. Numerical Implementation 

Details of our numerical finite-difference implementation 
can be found in m, but we review some of the key features 
here. 

We adopt a grid in three spatial dimension, using 
(N r , Nq, N<p) grid points. The grid is cell-centered, so that 
no grid points reside at r = 0 or sin 9 = 0. We use fourth- 
order differencing to evaluate most spatial derivatives (advec¬ 
tive terms are differenced with a third-order upwind scheme); 
this means that we need to pad the numerical grid with two 
layers of ghost zones. Except at the outer boundaries, where 
we impose simple outgoing-wave fall-off conditions, these 
ghost zones correspond to another zone in the interior grid 
(see Fig. 1 in da for an illustration). A ghost zone with coor¬ 
dinates 0 g and p g and a negative radius r g = — Ar/ 2 , for ex¬ 
ample, where Ar is the radial grid spacing, corresponds to the 
interior zone at 0 = 7r—0 g , (p = p g +it and r = —r g = Ar/2. 
The ghost zones can therefore be filled by copying function 
values from the corresponding interior zones. For tensor com¬ 
ponents, appropriate parity conditions have to be taken into 
account, since unit vectors in the ghost zone may point into 
the opposite direction of those in the corresponding interior 
zone (see Table I in |J5|). 

We implement a second-order version of the PIRK method 
for the time evolution. The stability of this method requires 
that the time step be limited by a Courant condition of the 
form 

At < CA min , ( 22 ) 

where C is a Courant factor and A m i n is the smallest distance 
between neighboring grid points. We evaluate this condition 
using simple coordinate distances, and chose C = 0.2 for 
all simulations in this paper. It is a well-known disadvantage 
of spherical polar coordinates that the accumulation of grid 
points in the vicinity of the origin leads to a severe limit on the 
time step. Nevertheless, we have performed all results shown 
in this paper with a serial implementation using uniform grids. 

In fl5l we have presented several tests of our code, includ¬ 
ing convergence tests for Teukolsky waves and single black 
holes. Because different parts of the code are differenced to 
different order, the order of convergence depends on which 
term dominates the error for the variable under consideration. 
In l27l we also used this code to simulate the collapse of non¬ 
linear gravitational waves to black holes. 

While our code does not make any symmetry assumptions, 
all simulations in this paper are axi-symmetric. We therefore 
choose the smallest possible number of grid points in the ip- 
direction, Ng> = 2, in all simulations presented here (but we 
refer to [15] for genuinely three-dimensional simulations with 
% > 2 ). 


B. Relativistic Hydrodynamics 

We have previously discussed an implementation of rela¬ 
tivistic hydrodynamics in spherical polar coordinates in ED. 
We briefly review our approach here, and also discuss new 
features of the approach used in this paper. 


1. A Reference-Metric Formulation of Relativistic Hydrodynamics 

The equations of relativistic hydrodynamics are based on 
the conservation of rest-mass, expressed by the continuity 
equation 

V a (p 0 u a ) = 0, (23) 

and conservation of energy-momentum 

V b T ab = 0. (24) 

Here V denotes the (four-dimensional) covariant derivative 
associated with the spacetime metric g^, po the rest-mass 
density, u a the fluid four-velocity, and T ab is the stress-energy 
tensor 

T ab = p 0 hu a u b +pg ab , (25) 

where h = 1 + e + p/pa is the enthalpy, p the pressure, and 
where e is the specific internal energy. We also define the 
Forentz factor between the fluid and normal observers as 

W = —n a u a = cm*. (26) 

The quantities po, p, e, and the fluid velocity v % , defined as 

v a = 2-ry a b u b = (0 ,u l /W + P l /a ) , (27) 

form the so-called primitive variables Q 

In many recent applications the above equations are brought 
into a flux-conservative form, so that high-resolution shock¬ 
capturing (HRSC) schemes can be used to find accurate nu¬ 
merical solutions. In the process, a new set of hydrodynami- 
cal variables, namely the conserved variables, are introduced. 
A particularly common such formulation is the so-called “Va¬ 
lencia” form l28l (see also (29, 30] for reviews). 

In the Valencia formulation, the continuity equation takes 
the form 

d t (e 6 ^D)+d j (f D ) j = 0, (28) 

the Euler equation is 

d t (e 6 ^Si) + dj{f s )i j = ( s s )i, (29) 


1 We note a typo in eq. (20) of ED, which should be replaced with eq. {27^ 

above. 



5 


and the energy equation becomes 

+ dj(f T y = s T . (30) 

Here 

D = W po (31a) 

Si = W 2 * hp 0 Vi (31b) 

T = W 2 hpQ—p~D (31c) 

are the density, momentum density and the internal energy as 
seen by a normal observer, 

(. f D ) j = ae 64, ^D{v j - /3 j /a) (32a) 

= ae 6<t> ^(W 2 hp 0 hv i (v j - ?/a) +pS i i \32b) 
(f T ) j = ae 6 ^ — ft /a) +pv-’) (32c) 

are the corresponding fluxes, and the two source terms are 

0»s)i = ae^^(-T 00 adia + T° k dip k + 

t jk d ajk / 2) (33a) 

s T = ae^yfitfiKij - (T 00 /3* + T oi )^a), (33b) 
where we have abbreviated 

t ij _ jiOOpipj + 2T 0i yi + T ij . (34) 


< |32| >. Alternatively one could, of course, work around this 
problem by simply scaling out these geometrical factors, but 
the reference-metric approach has other appealing features as 
well. All variables are now defined as tensors of weight zero 
(unlike those in the original Valencia formulation), and the 
formalism meshes well with the very similar approach used 
for the BSSN equations. The covariant derivatives of the spa¬ 
tial metric Vifjk = 4 e 4 < ^ 7 jkdip + e^Vijjk, for example, 


can be evaluated in terms of the derivatives V^jk that have 
already been computed for the connection coefficients ([7]). 


We evaluate the covariant derivatives V in the fluid equa¬ 
tions by expanding them in terms of connection coefficients; 
the continuity equation ([35]), for example, becomes 


+ dj(f D y = ~(f D y t%. (36) 


The appearance of the new terms on the right-hand side are 
a disadvantage of this approach. The vanishing of the right- 
hand side of the continuity equation in the original Valencia- 
formulation meant that, in an HRSC implementation, the total 
rest mass 


M 0 


J cfix^yau 1 po = 



(37) 


In curvilinear coordinates, the appearance of the determi¬ 
nant 7 in the above equations - in particular in the Euler equa¬ 
tion - poses a problem. Even for flat space in spherical polar 
coordinates we have 7 1 / 2 = r 2 sin$. The ^-dependence of 
this term leads to the appearance of non-zero terms on both 
sides of the Euler equation ( [29] ), even in spherical symmetry. 
Analytically these two terms cancel each other, but since in 
HRSC schemes the flux terms on the left-hand side are evalu¬ 
ated differently from the source terms on the right-hand side, 
numerical error will prevent a perfect cancellation. The re¬ 
sulting error breaks spherical symmetry, and can build up very 
quickly (see m for a more detailed discussion). 

In I2T1 we suggested a reference-metric approach, analo¬ 
gous to that applied to the BSSN equation in Section |Il A 1| 
to solve this problem]^] As we derived in ED, the resulting 
equations are very similar to those of the original Valencia 
formulation above, except that all appearances of ^/7 have 
to replaced with ^ 7/7 (which immediately solves the prob¬ 
lem discussed above), and all partial derivatives d have to be 
replaced with covariant derivatives with respect to the refer¬ 
ence metric, V. The continuity equation (28), for example, 
becomes 


dt (e 6 * ^hD) + ( f D y = 0. (35) 

From a computational perspective, the most important advan¬ 
tage of this reference-metric approach is that the geometri¬ 
cal factors r 2 sin 0 (and similar for other curvilinear coordi¬ 
nate systems) are eliminated from the definition of the fluxes 


2 We note that this problem has been recognized before. In general relativis¬ 

tic hydrodynamics this issue has also been addressed by 13111321 


is conserved exactly ; this is no longer the case in the reference- 
metric formulation (see Section [III B 2| below for a numerical 
example.) 

In f 2 H we therefore considered two different implemen¬ 
tations, a full approach, in which we adopted all three hy¬ 
drodynamic equations in the reference-metric version, and a 
partial approach, in which we adopted the Euler equation in 
the reference-metric version, but left both the continuity and 
the energy equation in the original Valencia formulation. In 
ED we found that the advantages of the partial approach out¬ 
weighed those of the full approach; however, our implemen¬ 
tation there did not adopt a factoring of tensor components. 


2. Factoring of Tensor Components 

In this paper we apply to all hydrodynamical variables the 
same factoring of tensor components that we described in Sec- 
tion |II A 2] The momentum density Si, for example, is written 
as 


Si=( rse \ , (38) 

V r sin 6s^ J 

(and similar for all tensorial variables), and the s r , sq and s^ 
are then evolved as the dynamical variables. We found that 
with this rescaling, and using the full approach, i.e. adopting 
the reference-metric approach for all hydrodynamical vari¬ 
ables, even a shock wave passing through the origin of the 
coordinate system will not lead to the formation of spikes or 
other numerical artifacts - we will present examples in Sec- 
tion lfflBl 
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3. Equation of State 


For the simulations presented in Section III B 


we construct 


the initial data using a poly tropic equation of state (EOS) 

P = npo, (39) 


where the polytropic constant k, is a measure of the entropy, 
and where T is related to the polytropic index n by r = 1 + 
1/n. The specific internal energy density e can then be found 
from the first law of thermodynamics to be 


1 P 


(40) 


During the dynamical evolution of these initial data we 
adopt a T-law EOS, meaning that we use eq. ( [40] ) to find P 
in terms of the conserved variables po and e. 

For all simulations in Section HUB I we will use units in 
which the polytropic constant n in (39) is unity, k = 1. How¬ 
ever, all dimensional quantities can be rescaled for any other 
value of k by recognizing that, in geometrized units, K n ' 2 has 
units of length. Density, for example, has units of inverse 
length squared, again in geometrized units. The star consid¬ 
ered in Section [Hi B 2| has a central density of p c = 0.2 in our 
units with n = 1 ; for any other value of k, the central density 
is then K~ n p c . 


4. Numerical Implementation 


We use a HRSC scheme to solve the equations of relativis¬ 
tic hydrodynamics in the above form (see, e.g. (33] 34’] for 
text-book treatments). In particular, we have implemented a 
second-order slope limiter reconstruction scheme, namely the 
monotonic centered limiter E3, to obtain the left and right 
states of the primitive hydrodynamic variables at each cell 
interface. We have also adopted the Harten-Lax-van-Leer- 
Einfeld approximate Riemann solver [ 36][37]]. We also refer to 
126, 381 for similar treatments and more detailed discussion. 

Following common practice we introduce an artificial at¬ 
mosphere to deal with the vacuum regions of spacetime, 
which would otherwise create numerical problems. We fol¬ 
low the prescription of (26l we set the density po to p atm = 
/ atm max(p 0 ) wherever it would otherwise be less than a 
threshold density p th res = /thresPatm, or wherever e ^ ^atm* 
Here we compute e at m from p atm using a polytropic EOS. 
For the simulations presented in Section|IIIB|we used / atm = 
10 -8 and /th res = 10 . 


III. OFF-CENTER SIMULATIONS 


t = 10.0M 



FIG. 1. The lapse function a in a slice of constant azimuthal angle 
<p, evolved on a grid of size (192,36, 2). The color-coded surface 
shows the data at time t m 10M, while the wireframe shows the 
initial data at t — 0 - the two surfaces are hardly distinguishable in 
the figure. The contour lines are drawn for a = j x 0.1, where j is 
an integer. 
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r/M 


FIG. 2. The lapse function a (top panel) and the conformal exponent 
<f> (bottom panel) as a function of radius r in the 6 = 0 direction, 
at t — 9.5M, for the same simulation as in Fig. [T] The insets show 
the difference of each function between t — 0 and t = 9.5M. The 
center of the black hole is located at r m 1 M. 


A. Vacuum 

1. Off- Centered Schwarzs child Black Holes 

As a first example of off-centered simulations we consider 
a Schwarzschild black hole that is placed away from the ori¬ 
gin of the coordinate system. We express this black hole in 


maximally-sliced “trumpet” coordinates (see (39)), which can 
be expressed analytically in parametric form (40) . In these co¬ 
ordinates, the black-hole horizon is located at a coordinate dis¬ 
tance rhor = 0.779M away from the black hole’s center. We 
evolve these data with the non-advective 1 +log slicing (17] ) 
and Gamma driver condition ([20]), using the analytical values 
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FIG. 3. In the top panel we show the value of the lapse a at the 
origin of the coordinate system, r — 0, as a function of time for 
a Schwarzschild black hole whose center is located at r — 1 M 
and 6 — 0. We show results for three different grid resolutions 
(7V96,1V18,2) with N = 1, 4/3 and 2. In the bottom panel we 
show the numerical errors A a rescaled with factors of N 2 , demon¬ 
strating at least second-order convergence. 


for the lapse and shift as initial data. Analytically, the result¬ 
ing evolution should result in all metric components remain¬ 
ing time-independent; any evolution away from the initial data 
is therefore caused by numerical error. 

We place the center of the black hole at r center = 1 M and 
0 = 0 , i.e. at a distance of 1 M from the origin of the coordi¬ 
nate system in the positive ^-direction (the direction towards 
the North Pole), and choose the numerical grid to extend to an 
outer boundary imposed at r max = 16M. We also choose the 
grid resolution such that the radial resolution A r = r max /7V r 
is similar to the angular resolution r A6 = vir/Ng at the center 
t center of the black hole, which implies 

Ng ~ N r . (41) 

T max 


sets is very small, so that they cannot be distinguished in the 
figure. In particular, we do not observe any problems arising 
at the origin of the coordinate system at r = 0. In Fig. [2] we 
show both the conformal exponent <p and the lapse function 
a at t = 9.5M along the axis pointing from the origin to the 
north pole, i.e. along 0 = 0. The insets in Fig. [2] display the 
difference between the values of these functions at t = 0 and 
t = 9.5 M. 

In the top panel of Fig. [3] we show the lapse function a , in¬ 
terpolated to the origin at r = 0 , as a function of time, for the 
three different resolutions N = 1, 4/3 and 2 . Analytically, the 
lapse should remain exactly constant; any departure from this 
constant value is therefore a measure of the numerical error. In 
the bottom panel we plot this numerical error Aa = a — a ana 
and multiply the result with N 2 . This graph demonstrates that 
the errors decrease slightly faster than second-order. 

Before closing this section we point out that, by placing the 
black hole away from the origin of the coordinate system, the 
numerical error becomes asymmetric about the center of the 
black hole. The angular resolution, r A0, at a given distance 
away from the center of the black hole, is smaller on the side 
of the black hole facing the origin of the coordinate system 
than on the side facing away from the origin (see Fig. 0 ) We 
have observed that this asymmetric error results in a slow drift 
of the black hole towards the origin. This drift cannot yet be 
seen in Figs. [T] and [2| but the fact that the lapse starts to de¬ 
crease at later times in Fig. [3j as the center of the black hole 
slowly approaches the origin, is a symptom of that drift. How¬ 
ever, this figure also demonstrates that this drifts converges 
away as the numerical resolution is increased. 


2. Head-on collision of two black holes 

As a second example of off-centered simulations in vacuum 
spacetimes we consider the head-on collision of two black 
holes from rest. As initial data we adopt Brill Lindquist m 
data, for which the spatial metric is conformally flat, the ex¬ 
trinsic curvature vanishes, and the conformal factor ip = is 
given by 




Mi M. 2 
r x r 2 


(42) 


In the following we present results for grids of size 
(7V96, M 8 , 2 ) with N = 1 , 4/3 and 2 . 

In Fig. [T] we show a “slice” of the lapse function a at a 
time t = 10M. Here and in the following, a slice, similar 
to a slice in an orange, represents the data as a function of r 
and 0 for a given value of the azimuthal angle p. Given that 
initial data for our simulations in this paper are axisymmetric, 
the particular value of p does not matter. We then graph the 
data as functions of Cartesian coordinates z = rcos# and 
x = r sin#. 

In Fig. [I] we show the data for our highest-resolution run 
with N = 2 (i.e. on a grid of size (192, 36, 2 )) at time t = 
10 M as the shaded surface. We also include the initial data at 
t = 0 as a wireframe representing each grid point used in this 
simulation. As expected, the difference between the two data 


Here r a = \x l — C l a \ is the coordinate distance to the center 
of the black hole located at C l a . We choose both black holes 
to start out on the axis. In order to obtain asymmetric data we 
choose M 2 = 2A4i, so that the total ADM mass M of the 
initial data is M = 3Mi. We place Mi at z\ = 3M\ = 
M, and M 2 at z 2 = — 6A4i = —2M. The center of mass 
is then at zqm = —M. After the two black holes coalesce, 
we expect them to merge close to their center of mass (we 
do not compute the energy or linear momentum emitted in 
gravitational radiation in these simulations.) This means that 
Mi will pass through the origin of the coordinate system prior 
to merger, making this a strong test of our numerical methods 
in spherical polar coordinates. 

We evolve these data with the advective version of 1+log 
slicing, eq. (18]), as well as the shift condition ([2T). As initial 
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t = 9.50M 
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FIG. 4. Snapshots for the head-on collision of two black holes. The 
top panel shows the initial data for the lapse function a , based on 
Brill-Lindquist data; the middle panel at the instant when the smaller 
black hole moves through the origin of the coordinate system; the 
bottom panel at a late time after merger. As in Fig. [T] the contours 
are drawn for values a = j x 0.1, where j is an integer. In these 
plots we show only every second grid point for clarity. 


data for these gauge conditions we use the “pre-collapsed” 
lapse 


a = ip ~ 2 = e -2 ^, (43) 

as well as vanishing shift, /3 Z = 0. We also experimented with 
an advective version of the Gamma-driver condition ( [20} , but 
these simulations crashed shortly after Mi passed through the 



FIG. 5. The lapse function a (top panel) and the conformal exponent 
4> (bottom panel), evaluated at the origin r — 0 of the coordinate 
system, as a function of time, for the head-on collision of two black 
holes. 


origin of the coordinate system. 

For the results presented in this paper we chose a numerical 
grid extending to an outer boundary at r max = 24Mi = 8M. 
Using expression © we chose the angular resolution simi¬ 
lar to the radial resolution at the initial position n = 3Mi 
of Mi. We present results for the three different grid sizes 
(N 128,7V48,2) with N = 1 , 3/2 and 2. 

In Fig. [4] we show slices of the lapse function a at three dif¬ 
ferent times. The top panel shows the initial data with the two 
black holes at their initial positions. The middle panel shows 
the data at a time of 9.5M, at the moment that M\ passes 
through the origin of the coordinate system. Finally, the bot¬ 
tom panel shows the remnant at a time well after merger. The 
grid lines in these slices represent every second grid point used 
in the simulations, and we note that the coordinate singular¬ 
ities at the origin of the coordinate system as well as on the 
axis do not cause any numerical problems. 

In Fig. [5] we show the lapse function a (top panel) and the 
conformal exponent 0 (bottom panel) interpolated to the ori¬ 
gin of the coordinate system as a function of time. The lapse 
function initially increases (until t ~ 4M), which is a result 
of the coordinate transformation from wormhole-type initial 
data to a trumpet geometry. It then drops to a value close to 
zero at t ~ 9.5 M, when M\ passes through the origin; the 
conformal exponent (j) has a sharp maximum at a similar time. 
At later times both the lapse and the conformal factor settle 
down to equilibrium values. 

We include results for the three different resolutions N = 1, 
3/2 and 2 in Fig. [5] These resolutions are sufficient to establish 
convergence at early times, before the origin of the coordinate 
system has become affected by the finite differencing across 
the singular conformal exponent at the center of black hole 
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FIG. 6. The rest-mass density po, the pressure p , and the specific 
internal energy e, interpolated to the axis, for the shock-tube problem 
described in the text. The solid line denotes the analytical solution 
at time t — 0.3; the crosses mark our numerical solution. The two 
different colors (blue and red) represent data on the northern and 
southern axis (i.e. for 0 = 0 and 0 — 7r); we do not observe any 
problem at the origin of the coordinate system, where the two colors 
meet. 


M\. At later times, however, these resolutions are not yet 
sufficient to establish convergence, as the error still appears to 
be dominated by higher-order error terms. We found a very 
similar behavior even at early times for only slightly smaller 
grid resolutions. 


B. Relativistic Hydrodynamics 


1 . Planar Shock-Tube 


t = 0.30 



s 


t = 0.30 



t = 0.30 



As a test of our implementation of relativistic hydrodynam¬ 
ics in spherical polar coordinates we first present results for a 
planar, special relativistic shock-tube problem. We consider a 
fluid at rest with two different homogeneous densities in the 
two hemispheres, separated by a diaphragm in the equatorial 
plane (i.e. at z = 0) until a time t = 0. At t = 0 the partition 
is removed, which results in a shock that propagates into the 
low-density region, while a rarefaction wave propagates into 
the high-density region. The analytical solution for this spe¬ 
cial relativistic shock tube problem is given in 1421 (see E3 
for the general solution). 


We choose an equation of state as described in Section 


IIB 3 in particular we set up polytropic initial data with a 


poly tropic index T = 1.2, and evolve these data with the 
Gamma-law EOS p0|). For the example presented here we 


FIG. 7. The rest-mass density po, the pressure p, and the specific 
internal energy e on a slice of constant azimuthal angle ip for a planar 
shock-tube problem. The color-coded surface shows the numerical 
data at time t — 0.3, while the wireframe shows the analytical data 
at the same time. The contour lines are drawn at eight equidistant 
values between the minimum and maximum values of the respective 
variables. 


set 


Po 


nr 8 , e < 7r/2 

io -7 , o > 7r/2 


(44) 


(in our units with k 1: see Section |II B 3[ >. While we do 
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evolve the spacetime together with the fluid for this test, we 
have chosen the densities sufficiently small so that the space- 
time remains very close to flat. 

For this test we imposed an outer boundary at r max = 0.5, 
and chose a numerical grid of size (192,96, 2). We use the 
“full” version of our reference-metric formulation of hydrody¬ 
namics in this Section, meaning that all hydrodynamical equa¬ 
tions are expressed in terms of a reference metric (see Section 


IIB 1). In the following we show results at time t = 0.3, after 
which the shock has travelled to z ~ 0.09. 


In Fig. [ 6 ] we show the analytical (solid line) and numerical 
solutions (crosses) at t = 0.3. For this graph we interpo¬ 
lated the numerical data to the axis 6 = 0 , so that they can be 
compared directly with the analytical solution. For most parts 
of the solution the agreement is excellent. As expected, the 
shock itself is spread out over about three grid points. Also 
as expected, the contact discontinuity (at about z = 0.06 in 
Fig. [ 6 ]) poses the greatest numerical challenge, especially in 
the specific internal energy. However, we have compared with 
lower-resolution results to confirm that increasing the numer¬ 
ical resolution results in an improved representation of this 
discontinuity. What distinguishes this figure from many other 
shock-tube tests, however, is the presence of the origin of the 
coordinate system. In order to highlight this, we plotted grid 
points on the northern axis (6 = 0 ) in red, and those on the 
southern axis (0 = i r) as blue. We do not observe any numeri¬ 
cal problems at the origin of the coordinate system, where the 
two colors meet. 


Fig. [7] shows the fluid variables at the same time t = 0.3, 
but represented as surface plots. The colored surfaces show 
numerical “slices” of the fluid variables, with the grid lines 
again representing all grid points used in these simulations. 
The (square) wireframe shows the analytical solution. We 
again notice very good agreement between the numerical and 
analytical solutions. The shock is smeared out across three 
grid points everywhere; however, since the angular resolution 
decreases further away from the axis (i.e. for larger values of 
x) 9 the shock also becomes less sharp. The contour lines in¬ 
dicate that the entire shock remains very close to planar, even 
though it is represented in a coordinate system that does not 
reflect this symmetry. Most importantly, we again observe that 
no numerical problems arise at the origin of the coordinate 
system. This test therefore demonstrates that our numerical 
implementation, which includes the reference-metric formu¬ 
lation of hydrodynamics together with a proper factoring of 
all tensor components, allows for the simulation of a shock 
wave in spherical polar coordinates, even at the origin of the 
coordinate system. 

We also experimented with planar shocks that originate 
from a partition located at values of z < 0 (i.e. not in the 
equatorial plane). The disadvantage of this setup is that the 
initial discontinuity is not aligned with cell boundaries of our 
numerical grid. This results in a noisy representation of the 
initial data, which, in turn, results in numerical errors that 
are larger than those encountered in the experiments described 
above - even before the shock wave reaches the origin. While 
we did not encounter any instabilities or problems at the ori¬ 
gin in these simulations, even after the shock wave had passed 



FIG. 8. The rest-mass density po on a slice of constant azimuthal an¬ 
gle p. The color-coded surface shows the data at time t = 175.4M, 
while the wireframe shows the initial data at t = 0. The contour 
lines are drawn for eight equidistant densities. 


through it, we decided to focus on shocks originating from 
the equatorial plane here, so that the numerical evolution is 
not affected by numerical noise in the initial data. 


2. Off-centered Tolman-Oppenheimer-Volkoff Stars 


As a last example we consider a static equilibrium star, i.e. a 
solution to the Tolman-Oppenheimer-Volkoff (TOV) equa¬ 
tions |43j 44 f|. We set up the initial data as a relativistic poly¬ 


trope with T = 2 (see Section IIB 3). For this poly tropic 
index, the maximum-mass model has a central rest-mass den¬ 
sity of p^ ax = 0.318, a rest-mass of M 0 = 0.180, and an 
ADM mass of M = 0.164 (in our code units with n = 1). 
For the simulations shown in this Section we choose a star 
with central density of p™^ = 0.2, for which the rest mass 
is M 0 = 0.172 and the ADM mass is M = 0.157. The areal 
radius of this star is R = 0.866 = 5.52M, while the isotropic 
radius is r = 0.700 = 4.46M. We evolve this star with non- 
ad vective 1+log slicing ( fTT) and zero shift, so that, analyti¬ 
cally, the star should remain static. Any deviation from the 
initial data is therefore a measure of the numerical error. 

We place the center of this spherical star at r = 0.5 = 
3.18M on the northern axis (0 = 0), so that the origin of the 
coordinate system is inside the star (see Fig. [8]). We evolve 
this stellar model on a grid that extends to an outer boundary 
at r max = 3 = 19.1M. We use ([41) to match the angular res¬ 


olution to the radial solution at the center of the star; in the fol¬ 
lowing we show results for grid resolutions of (TV 64, TV32, 2 ) 
with TV = 1, 3/2 and 2. 

The colored surface in Fig. [ 8 ] shows a slice of the rest- 
mass density po at time t = 175.4M, evolved on the highest- 
resolution grid with TV = 2. The contour lines demonstrate 
that the star remains nearly spherical. Also included in the 
figure is a wireframe that shows the initial data. It can barely 
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FIG. 9. The central value of the rest-mass density po, interpolated to 
r iss 0.5 ft= 3.18M and 0 — 0 for a star with initial maximum density 
p™ ax = 0.2. The top panel shows the data themselves for different 
resolutions; the bottom panel shows the relative errors Ap 0 /p™ ax , 
with Ap 0 = po — p^ ax , rescaled with a factor N 2 . We observe 
second-order convergence until t ~ 25 M, at which point the center 
of the star is affected by errors originating from the surface of the star. 
The simulations converge even after that time, but at a rate slightly 
less than second-order. 


be noticed in the figure that the star has slightly drifted to¬ 
wards the origin - a numerical artifact very similar to that de¬ 
scribed in Section III A 1 for a single black hole. Again, this 
effect becomes smaller with increasing resolution. The ori¬ 
gin of the coordinate system can be seen in the envelope of 
the star; clearly, the presence of the coordinate singularities at 
the origin and on the axis do not cause any problems in this 
simulation. 

In Fig. [9] we show Po enter , the rest-mass density po interpo¬ 
lated to the center of the star at r = 0.5 = 3.18M and 6 = 0. 
Analytically, this value should remain at p™* = 0.2, but nu¬ 
merical error causes small deviations. In the top panel of Fig|9] 
we show po enter f° r the three different resolution TV = 1, 3/2 
and 2 as a function of time, in the bottom panel we show the 
relative errors (po enter — Po Lax )/p 1 o ax multiplied with factors 
TV 2 . Our results show second-order convergence until a time 
of t ~ 25M; after that, the center of the star is affected by 
numerical errors originating from the stellar surface, where 
discontinuous derivatives of fluid and metric variables lead to 
a slower rate of convergence. 

In Fig. [TO] we show similar results for the total rest mass 
Mo, given by the integral ( [37] ). In the top panel we show nu¬ 
merical values of M 0 as a function of time for the three differ¬ 
ent resolutions TV = 1, 3/2 and 2; in the bottom panel we show 
the relative numerical errors multiplied with TV 3 . These find¬ 
ings suggest approximately third-order convergence of the rest 
mass. As discussed in Sect. |nm the origin of this error is 
the new term on the right-hand side of eq. ([36]), which appears 


FIG. 10. The rest mass Mo as a function of time for the same sim¬ 
ulations as Fig. [9] The top panel shows the data for different resolu¬ 
tions; the bottom panel shows the relative errors AMo/M 0 ana , with 
AMo = Mo — Mo na , rescaled with a factor TV 3 . We find that the 
rest mass converges to approximately third order. 


in the “full” reference-metric formulation of hydrodynamics. 
In Fig. [IT] we therefore compare results for both the rest-mass 
density and total rest-mass as obtained in the “full” and the 
“partial” approaches. The top panel of Fig. El shows that 
the numerical errors for the rest-mass density are somewhat 
smaller in the full approach than in the partial approach, while 
the bottom panel demonstrates that, as expected, the total rest 
mass is conserved much better in the partial approach (the 
only error originates from our treatment of the atmosphere, 
see Section |nB4| ). 


IV. SUMMARY AND DISCUSSION 


We have recently developed a new approach for numerical 
relativity simulations in spherical polar coordinates (T5j 2 J ]. 
Our approach is not based on a regularization of the equa¬ 
tions, and instead deals with the coordinate singularities at the 
origin and on the axis of the coordinate system by using a 
reference-metric formulation of the equations, a factoring of 
all tensor components, and a partially implicit Runge-Kutta 
method. Unlike previous approaches that employed a regu¬ 
larization of the equations, our method does not rely on any 
symmetry assumptions. 

While we have previously presented several tests of these 
methods, most of these test cases featured a symmetry about 
the origin. In this paper we therefore present new “off- 
centered” simulations that highlight the stability of our meth¬ 
ods in the absence of such a symmetry. We perform several 
stringent tests to show that the coordinate singularity at the 
origin of the coordinate system does not pose any computa- 
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FIG. 11. A comparison between the full and partial approaches, for 
the resolution N = 3/2 shown in Figs. [5] and 10 The top panel 
shows that the errors in the central density are somewhat smaller in 
the full approach than in the partial approach. The bottom panel 
shows that, as expected, the rest mass in the partial approach is con¬ 
served almost exactly; the only error in the rest-mass arises from the 
treatment of the atmosphere (see Sects. inB~n an d [HB4l ) 


tional problems - even for a black hole that drifts through the 
origin. 

The other purpose of this paper is to discuss an alternative 
implementation of the reference-metric formulation of rela¬ 
tivistic hydrodynamics. Unlike in [2T), we here apply the 


same factoring of hydrodynamical tensor components as we 
use for Einstein’s field equations. We perform several simula¬ 
tions to test and calibrate this implementation. Perhaps most 
importantly, we demonstrate that our method is able to fol¬ 
low the propagation of a shock wave through the origin of the 
spherical polar coordinate system. We are not aware of a pre¬ 
vious solution to this problem. 

We also discuss the respective advantages of the “full” and 
“partial” approach in our reference-metric formulation of rela¬ 
tivistic hydrodynamics (see Sect. IIB 1| ). The partial approach 
applies this formulation only to the Euler equation, while the 
full approach applies it to all hydrodynamical variables and 
equations. The advantage of the partial approach is that the 
right-hand side of the continuity equation vanishes, so that, 
in HRSC implementations, the rest mass is conserved exactly 
(except for errors originating from the treatment of the atmo¬ 
sphere). The full approach, on the other hand, allows for a 
more accurate treatment of the origin of the coordinate sys¬ 
tem, and avoids any spiky or singular behavior even for a 
shock-wave passing through the origin. Which one of the two 
approaches is preferable may therefore depend on the appli¬ 
cation. 
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